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Multiparticle collision dynamics (MFC), a particle-based mesoscale simulation technique for com¬ 
plex fluid, is widely employed in non-equilibrium simulations of soft matter systems. To maintain a 
defined thermodynamic state, thermalization of the fluid is often required for certain MFC variants. 
We investigate the influence of three thermostats on the non-equilibrium properties of a MFC fluid 
under shear or in Foiseuille flow. In all cases, the local velocities are scaled by a factor, which is 
either determined via a local simple scaling approach (LSS), a Monte Carlo-like procedure (MCS), 
or by the Maxwell-Boltzmann distribution of kinetic energy (MBS). We find that the various scal¬ 
ing schemes leave the flow profile unchanged and maintain the local temperature well. The fluid 
viscosities extracted from the various simulations are in close agreement. Moreover, the numerically 
determined viscosities are in remarkably good agreement with the respective theoretically predicted 
values. At equilibrium, the calculation of the dynamic structure factor reveals that the MBS method 
closely resembles an isothermal ensemble, whereas the MCS procedure exhibits signatures of an adi¬ 
abatic system at larger collision-time steps. Since the velocity distribution of the LSS approach is 
non-Gaussian, we recommend to apply the MBS thermostat, which has been shown to produce the 
correct velocity distribution even under non-equilibrium conditions. 


I. INTRODUCTION 

Multiparticle collision dynamics (MPG) is a particle- 
based mesoscale simulation technique for fluids which has 
been introduced about a decade ago m- By now, it has 
been developed into one of the major simulation tech¬ 
niques for complex fluids and has been applied to a broad 
range of soft matter systems. Examples cover equilibrium 
colloid [DIHE] and polymer [aiiiiiiiHTi] solutions and, 
more importantly, non-equilibrium systems such as col¬ 
loids [THlfeS] , polymers [T^ , vesicles [SS] , and cells 

[37l|38j in flow fields, colloids in viscoelastic fluids [39H4T] . 
as well as of self-propelled spheres , rods [31 US] , 

and other microswimmers |46H50] . Moreover, extensions 
have been proposed to fluids with non-ideal equations of 
state m and mixtures [52]. 

The MFC algorithm consists of two discrete steps— 
streaming and collision—, and shares many features with 
the Direct Simulation Monte Carlo (DSMC) approach 
[53| . Although space is discretized into cells to define 
the multiparticle collision environment, both the spa¬ 
tial coordinates and the velocities of the particles are 
continuous variables. Hence, the algorithm exhibits un¬ 
conditional numerical stability and satisfies a H-theorem 
[HEiiMi. 

MFC refers to a class of algorithms, which differ in 
their collision rules [4]. In the original version of MFC, 
denoted as stochastic rotation dynamics (SRD) [T1155]. 
collisions consist of a stochastic rotation of the relative 
velocities of the particles in a collision cell. Other al¬ 
gorithms assign Maxwellian distributed random relative 
velocities to the particles in a collision cell at every colli¬ 
sion step, such that the momentum of the collision cell is 
conserved (MFC-AT) [8l[T9l[56]. MFC defines a discrete¬ 
time dynamics which has been shown to yield the correct 
long-time hydrodynamic behavior [H [^. A conse¬ 


quence of the discrete dynamics is that the transport co¬ 
efficients depend explicitly on the collision-time interval 
[m[ssi[iMi], which in turn permits control of fluid 
properties. 

In many non-equilibrium systems, temperature has to 
be controlled to ensure a stationary state. A defined tem¬ 
perature is inherent in the MFC-AT approach [idiins], 
but an additional mechanism has to be provided for the 
MFC-SRD version, since it conserves energy. Various 
constant temperature simulation schemes have been pro¬ 
posed [33H75] : not all of them ensure that a canonical 
ensemble is achieved and not all of them conserve mo¬ 
mentum. 

Under equilibrium conditions, momentum can be con¬ 
served by velocity scaling schemes [HEKMIEIKTSI. In its 
simplest form, velocity scaling keeps the kinetic energy 
of a system at the desired value by multiplying the veloc¬ 
ities of all particles by the same factor m- This corre¬ 
sponds to an isokinetic rather that an isothermal ensem¬ 
ble. As far as MFC is concerned, a local cell-level scaling 
scheme (LSS) can be implemented, where a scale factor 
is determined for every cell independently. To achieve a 
canonical distribution of kinetic energies, more sophisti¬ 
cated cell-level approaches have been proposed based on 
a Monte Carlo criterion [aEniTg, which we denote as 
Monte Carlo scaling (MCS), or by exploiting the appro¬ 
priate distribution of the kinetic energy (Gamma distri¬ 
bution), denoted as Maxwell-Boltzmann scaling (MBS) 

m- 

As is well known, under non-equilibrium conditions an 
inappropriate thermostat may introduce a bias into sys¬ 
tems with an a priori unknown velocity profile |77] . To 
prevent a bias in MFC simulations, the relative parti¬ 
cle velocities with respect to the center-of-mass velocity 
of a collision cell are scaled, which yields a profile un¬ 
biased thermostat (FUT) [501 [7S] and renders a global 
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scaling scheme inappropriate. However, recent detailed 
MFC simulation studies with a particular collision rule 
indicate a substantial interference of certain thermostats 
with the flow field [75]. The comparison of viscosities, ex¬ 
tracted from the parabolic flow profiles of Poiseuille flows, 
yields surprisingly large deviations between the values 
extracted from simulations and those determined theo¬ 
retically. To avoid such an effect. Ref. [75] suggests to 
exclude the flow and shear directions from thermostat¬ 
ting. Since, the large deviations are rather unexpected, 
we perform non-equilibrium MPC-SRD simulations with 
the “standard” three-dimensional collision rule in order 
to unravel the underlying cause. 


A suitable thermostat is essential for accurate and re¬ 
liable simulation results, and thermostats failing for flow 
fields like shear or Poiseuille flow are obviously inade¬ 
quate for more complex flows. Hence, studies of the reli¬ 
ability of a thermostat combined with a particular MPC 
collision rule are important. In this article, we character¬ 
ize the influence of the LSS, MCS, and MBS thermostats 
on the velocity profile of simple shear and Poiseuille flow. 
Thereby, we determine the viscosity for various collision 
times and compare it with the theoretical expression. As 
a reference, we also determine the fluid viscosity at equi¬ 
librium via the transverse fluid-velocity autocorrelation 
function in Fourier space without any thermostat m- 
We find very good agreement between the velocity pro¬ 
files for the various thermostats. The differences between 
the relative viscosities are below 2% and thus agree with 
each other within the accuracy of the simulations. The 
analytically derived expression of the viscosity of a MPC 
fluid is based on some approximations, specifically the 
molecular-chaos assumption. Hence, it is not a priory ev¬ 
ident to which extent it describes simulation results. We 
find a remarkably good agreement between the viscosity 
extracted from simulations and the theoretical expres¬ 
sion, which emphasizes the importance of the theoretical 
result. Moreover, we calculate the equilibrium dynamic 
structure factor of fluids thermalized by the various ther¬ 
mostats. For systems with temperature controlled by 
LSS or MBS, the structure factor lacks a Rayleigh line, 
which corresponds to an isothermal ensemble [7i [HOj. 
For systems with a MCS thermostat, a Rayleigh peak 
appears for larger collision-time steps, hence, no isother¬ 
mal system is simulated. 


The article is organized as follows. In Sec.|TT]the model 
and simulation approach are introduced. Section HI 


describes the various thermostats. Results for equilib¬ 
rium hydrodynamic fluctuations, specifically transverse 
velocity autocorrelation functions and dynamics struc¬ 
ture factors, are presented in Sec.|T^ The results of non¬ 
equilibrium simulations are presented in Sec. jV] and vari¬ 
ous implications are discussed in Sec. |VI| Finally, Sec. m 
summarizes our findings. 


II. MODEL AND METHODS 
A. Multiparticle Collision Dynamics 

We consider a MPC fluid of N point particles of mass 
m. Without external field, the particles move ballisti- 
cally during the streaming step, i.e., their positions 
are updated according to 

ri{t + h) = r^{t)+ hvi{t), ( 1 ) 

where Vi{t) is the velocity of particle i at time t and 
h is the collision-time step. In the collision step, the 
simulation box is partitioned into cubic collision cells of 
side length a. In the SRD version of MPC, the relative 
velocity of each particle, with respect to the center-of- 
mass velocity of the particular cell, is rotated by a fixed 
angle a around a randomly oriented axis. Hence, the 
velocity after a MPC step is 

Vi{t + h)= Viit) + [R(a) - E] {vi{t) - Vcm{t)), (2) 

where R(a) is the rotation matrix, E is the unit matrix, 
and 


1 , ^ 

'^j I ( 3 ) 

i=i 

is the center-of-mass velocity of the Nc particles con¬ 
tained in the cell of particle i mi- The random ori¬ 
entation of the rotation axis is chosen independently at 
every collision step and for every cell |75] . To insure 
Galilean invariance, a random shift is performed at every 
collision step m- In a collision step, mass, momentum, 
and energy are conserved which leads to the build-up of 
correlations in the particle motion and gives rise to hy¬ 
drodynamic interactions. 


B. Boundary conditions 

We typically apply three-dimensional periodic bound¬ 
ary conditions with a cubic simulation box of side length 
L and volume V = . In many systems, e.g., in simula¬ 

tions of Poiseuille flow, solid walls are present, commonly 
with no-slip boundary conditions. The discretization into 
collision cells requires particular measures to ensure the 
no-slip condition. We follow the approach suggested in 
Ref. [I8| . which applies the bounce-back rule, a random 
shift of the collision lattice in all spatial directions, and 
partial filling of surface cells by phantom particles. For 
systems with parallel walls, the random shift perpendic¬ 
ular to the walls is implemented as follows [5T]. Without 
random shift, collision-cell boarders are chosen to coin¬ 
cide with the respective wall. To enable a random shift, 
an additional layer of empty collision cells is added in 
one of the walls. In a random shift, the whole collision 
lattice is then shifted by a uniformly distributed random 
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displacement G [0,a]. The typically appearing partially 
occupied cells at the walls cause a violation of the no-slip 
boundary condition, since the average fluid velocity par¬ 
allel to the surfaces is no-longer zero in the surface cells 
m- To restore no-slip boundary conditions, usually a 
phantom particle is added to every cell intersected by a 
wall and occupied by fluid particles smaller than the 
average number of particles (iVc), such that the average 
particle density is restored. However, this does not com¬ 
pletely prevent slip, because the average center-of-mass 
position of all particles in a collision cell—including the 
phantom particle, which is placed in the center of the 
wall-occupied part of the collision cell—does not coincide 
with the wall. In order to fully account for the no-slip 
boundary condition, we adopt the following modification 
of the original approach m- To treat a wall cell on the 
same basis as a cell in the bulk, i.e., the number of fluid 
particles satisfies a Poisson distribution with the average 
{Nc), we take fluctuations in the particle number into 
account by adding Ngp particles to every cell partially 
occupied by a wall such that {Ngp + Ngc) = {Nc). The 
momentum P of all phantom particles of a cell is taken 
from the Maxwell-Boltzmann distribution with the vari¬ 
ance mNgpkBT and, at equilibrium, zero average. There 
are various ways to determine the number Ngp. For a 
system with two parallel walls, we suggest to use the 
number of fluid particles in the surface cell intersected 
by the opposite wall. The average of the two numbers 
is equal to (Nc). Alternatively, (Ngp) can be taken from 
a Poisson distribution with average (Nc) accounting for 
the fact that there are already Ngc particles in the cell. 
Collisions are then performed with all the particles in 
the cells. The center-of-mass velocity of the particles in 
a boundary cell is 

Instead of a single phantom-particle momentum and a 
single extra collision-lattice layer, an additional collision¬ 
cell layer can be added in every wall and explicitly 
be filled with randomly placed phantom particles with 
Maxwellian distributed velocities in every collision step. 
If the layers are sufficiently large, the fluctuations of the 
particle number in a collision cell are close to that of a 
bulk cell. 

Other approaches for no-slip boundary conditions have 
been suggested and analyzed [i El uni Hi isn- Some of 
them do not include phantom particles. However, based 
on our experience, we consider the approach with phan¬ 
tom particles as most appropriate for no-slip boundary 
conditions, which yields the correct hydrodynamic be¬ 
havior not only for solid walls, but also for solid particles 
immersed in a MFC fluid. 

The presence of external fields may require an adap¬ 
tation and modification of the boundary conditions to 
ensure the no-slip requirement. For simple shear such a 
adjustment is described in Ref. [ST], and for Poiseuille 


flow in Ref. EHj. We will come back to this aspect in 
Sec.lTml 

C. Shear Flow 

Shear flow is implemented by Lees-Edwards periodic 
boundary conditions [52J [51] , which yields a linear veloc¬ 
ity profile m- The shear viscosity ry of the fluid follows 
from the stress tensor via 77 = (Jxyli for shear along 
the cc-axis and the gradient direction along the j/-axis of 
the Cartesian reference frame; 7 denotes the shear rate. 
As shown in Ref. [61], the instantaneous stress of the 
MFC fluid is given by 

1 Af 7/1 ^ 1 ^ 

'^'^v = -yYl ~Vh'^ 

i—1 i—1 i—1 

(5) 

hence, a^y = {<^xy)- Via, ol G {x,y,z}, is the velocity 
after streaming, but before collision, whereas Via is the 
velocity after collision, but before streaming, and l^Pia 
is the momentum change of particle i during a collision. 
Here, Vi is the position of the particle in the grid-shifted 
frame. Note that due to Lees-Edwards boundary con¬ 
ditions, all particles are inside of the primary periodic 
box and the velocities along the flow direction are corre¬ 
spondingly adjusted [STj . 

D. Poiseuille Flow 

For the Poiseuille flow simulations, we confine the fluid 
between two solid walls, which are parallel to the xz- 
plane of the Cartesian reference frame, and apply peri¬ 
odic boundary conditions along the x- and z-axis. Flow 
is induced by a constant force mg acting on every fluid 
particle. Therefore, the particle velocities and positions 
are updated according to 

Vix{t-G h) = Vix + gh, ( 6 ) 

fix{t + h) = rix{t) + Vix{t)h+]^gh?' (7) 

along the flow direction. Note that the circumflex in¬ 
dicates quantities after streaming but before collision. 
To satisfy the no-slip boundary condition, we apply the 
bounce-back rule during streaming and take into account 
phantom particles during collisions. We consider two 
variants for the calculation of a phantom-particle mo¬ 
mentum. In the simpler version, the average momentum 
{P) is set to zero. However, this implies a residual slip. 
Following the proposition for shear flow in the presence 
of walls in Ref. m, we assign a finite mean velocity to 
every phantom particle according to its position in the 
collision cell relative to the wall-fluid interface. Thereby, 
we place a phantom particle in the center of the part of 
the collision cell inside a wall. The average velocity itself 
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is determined by the desired parabolic flow profile. A 
similar approach has been adopted in Ref. m- 

E. Viscosity 

Analytical expressions for the MFC fluid viscosity have 
been derived applying various methods [SHll |55j [58l - 
[SI1IH2]- In general, the viscosity r] = + r]'^ comprises a 

kinetic contribution due to streaming of the fluid par¬ 
ticles, and a collisional contribution j]^. For the latter, 
an exact expression can be derived, which is given by 


h = St are considered to cover the collisional-dominated 
as well as the streaming-dominated regime. The size of 
the cubic simulation box is set to L = 30a if not oth¬ 
erwise stated. For the shear flow simulation the choose 
the shear rate yr = 10“^, and for the Poiseuille flow 
simulations, we set g = 10“^a/T^. 

For an efficient simulation of the MFC fluid dynamics, 
we exploit a graphics processor unit (GFU) based version 
of the simulation code [HH| . 

III. THERMOSTAT 


1 


C 


Nma? 

18Vh 


(1 — cos(a)) 



( 8 ) 


for (Nc) S> 1 in three dimensions. Here and in the fol¬ 
lowing, we neglect fluctuations in the particle number 
in a collision cell, which is justified for average particle 
numbers (Nc) > 5, since we omit a term of order 
Due to correlations in the particle velocities, the kinetic 
contribution can only be derived within certain approxi¬ 
mations. Applying the molecular chaos assumption, i.e., 
velocity correlations between different particles are ne¬ 
glected, the kinetic contribution is 


V 


k 


NkeTh r 5 (A^c) _ 

2V ((Ac) — 1)(2 — cos(a) — cos(2a)) 

(9) 


Evidently, the collisional contribution dominates for 
small and the kinetic one for large collision-time steps, 
which corresponds to a fluid-like behavior in the first case 
and gas-like behavior in the second case as expressed by 
the Schmidt number |84| . 

In Ref. [84] and especially in Ref. [85| for two- 
dimensional systems, improved analytical expressions are 
provided for taking into account velocity correlations. 
It is important to note that fluid correlations yield a 
somewhat larger value than that predicted by the 
molecular-chaos assumption. 

The total viscosity ry = 77 ^ -|- is evidently dominated 
by for —)■ 00 and ry'^ for —>■ 0. Since g‘^ is calculated 
without any approximation, g provides an exact descrip¬ 
tion for ft, —0. Moreover, the applied molecular chaos 
assumption applies well for ft —>■ 00 . Hence, g is well 
described quantitatively by the theoretical expression in 
both limits. 


F. Parameters 

All simulation are performed with the rotation angle 
a = 130°, and the mean number of particles per collision 
cell {Nc} = 10. Length and time are measured in units 
of the collision cell size a and t = yjma^/{ksT), respec¬ 
tively, where T is the temperature and ks the Boltzmann 
constant. Various collision times between ft = O.lr and 


We perform simulations applying the thermostats 
mentioned in the introduction, namely local simple scal¬ 
ing (LSS) [ 73 , Monte Carlo scaling (MCS) as suggested 
in Ref. |B], and Maxwell-Boltzmann scaling (MBS) [75j . 
In all case, the relative velocities Avi = Vi — Vcm of the 
particles in a collision cell are scaled by a constant factor 
which typically differs for every cell and collision-time 
step. Hence, the relative velocities after collision Av{ 
are given by Av{ = ^Avi. Since the total relative 
momentum of a collision cell is zero, such a scaling leaves 
the total momentum of a cell unchanged. 

Simple Scaling —In the simple scaling approach, the 
scale factor ^ is chosen as 




l3{Nc-l)kBT 


2Ek 


with the kinetic energy of a collision cell 


Ek 


-| * c 


( 10 ) 


( 11 ) 


The factor Ac — 1 accounts for the fact that Ek is 
calculated in the center-of-mass reference frame of a col¬ 
lision cell. Strictly speaking, LSS conserves the average 
kinetic energy rather than the temperature |75j . This 
implies that the energy fluctuations are incompatible 
with that of an isothermal ensemble and the distribution 
of velocities in a collision cell is not Maxwellian. 


Monte Carlo Scaling —We implement the Monte Carlo 
scaling method in the version proposed in Ref. |5] , which 
satisfies detailed balance in contrast to earlier versions 
[65] . In brief, a factor e is randomly chosen in the inter¬ 
val [0,^] and ^ is set to either I -|- e or 1/(1 -|- e), each 
with the probability 1/2. The velocity scaling itself is 
performed following a Monte Carlo-type criterion, with 
the probability pA = min[l. A], where [HllT^ 

A = exp (-(^ - l)Ek/kBT) . (12) 

The choice of C S [0.05,0.3] and the frequency of scaling 
determine the relaxation time to approach the desired 
temperature T. We set ( = 0.1. The method has 
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been shown to yield the correct velocity distribution 
and has successfully been applied in various simulation 
studies PlIHlIHTHHn]- 


Maxwell-Boltzmann Scaling —In the Maxwell-Boltz- 
mann scaling approach, the known distribution of the 
kinetic energy of the MFC ideal-gas particles is exploited 
to determine the scale factor m- Thereby, the distribu¬ 
tion of the kinetic energy is given by (F distribution) 


Here, / = 3{Nc — 1) is the number of degrees of free¬ 
dom of the fluid particles in the considered collision cell, 
and r(a;) is the gamma function. In the limit / —>■ oo, 
the F distribution turns into a Gaussian function with 
mean (Ek) = fkBT/2 and variance f{kBT)^/2. At ev¬ 


ery collision, a random kinetic energy is taken from 
the distribution function (13) for every collision cell and 
the respective scale factor for the velocities is set to 


C = yj Ek/Ek, 


(14) 


with the kinetic energy E^ of Eq. (11). Taking the aver¬ 


age of the kinetic energy after scaling, we obtain 


T Na /.2 \ 

{Ek) = 2 H = ( y H 

i=l \ ^ / 


= (Ek 


(15) 


Hence, the average of the kinetic energy of a collision 
cell is equal to the desired mean of the distribution func¬ 
tion (13). Note that in Ref. [75], a different scale factor 


is provided, which may simply be a misprint, otherwise 
the factor would not provide the correct average kinetic 
energy. 

As has been shown in Ref. ITS], the MBS approach 
yields the correct distribution function of the particle ve¬ 
locities at the collision cell level, even for strong external 
fields, whereas LSS fails even at equilibrium. 

When phantom particles are taken into account, the 
MFC particles next to a wall are thermalized by the 
phantom particles. For sufficiently weak (external) fields, 
this energy exchange suffices to control the temperature 
in the whole system. For strong fields, the energy trans¬ 
port is not fast enough to ensure the desired temperature 
over the whole system m- Here, one of the additional 
thermostating schemes has to be applied. 



FIG. 1. (Color online) Transverse velocity autocorrelation 
functions of a MFC fluid for the collision-time steps hjr = 
0.1, 1.0, 0.2, and 0.5 (left to right). The fc-values are k = 
2'KnlL, with L = 30a and n = 1, 2, 3, and 4. The fit of the 
exponential function 0 (dashed lines) yields the kinematic 
viscosities presented in Table |I| 


correlation functions, on the one hand, to extract the 
fluid transport coefficients from equilibrium velocity au¬ 
tocorrelation functions, and, on the other hand, to verify 
the kind of simulated ensemble. 

Within the linearized Navier-Stokes equations 0110], 
the transverse hydrodynamic (shear) modes are indepen¬ 
dent of the longitudinal (acoustic, entropy) modes [75] . 
For an adiabatic system, i.e., a system in which the en¬ 
ergy of the fluid is locally conserved, the longitudinal 
modes are coupled. In particular, the temperature (or 
entropy) fluctuations are coupled to the density and lon¬ 
gitudinal velocity fluctuations O EO]. In contrast, in 
an isothermal system, temperature fluctuations are sup¬ 
pressed and controlled by the (local) thermostat, which 
implies a decoupling of the density and velocity fluctua¬ 
tions from the equation of the temperature fluctuations. 
The respective modifications of the transport properties 
are reflected in the density autocorrelation function, e.g., 
the dynamic structure factor S(k,io). 

In the following two subsections, we will address the 
transverse velocity correlation function and the density 
fluctuations via the dynamic structure factor. 


IV. RESULTS: HYDRODYNAMIC 
FLUCTUATIONS 

The hydrodynamic properties of a MFC fluid coincide 
with those of the linearized fluctuating hydrodynamic 
equations for sufficiently large length and time scales 
[D n EH ET] m Eg. Hence, we can use hydrodynamic 


A. Transverse Velocity Correlation Fnnction 

The MFC fluid viscosity can be determined by equilib¬ 
rium simulations, independent of any thermostat, from 
the transverse velocity correlation function in Fourier 
space [M1IS3ES]. For the considered periodic system, 
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the velocity in fc-space is defined as 


is given by m 


N 

(16) 

i=l 

with kjj = 27rn^/L, j3 £ {x,y,z'\, G Z, and k ^ 
0. From the linearized hydrodynamic equations [STlEi 
IHOIISIISS], the normalized transverse velocity correlation 
function Cy{k,t) is obtained, which decays exponential 
as 


27rS'(fc, to) = 


7 — 1 2DTk^ 

7 + {Dxk'^)^ 


( 20 ) 




. (w + Csk)^ + (r«fc 2)2 (cu - Cskf + (r,fc2)2 

+ ^ [r, + (7-i)i^T]- 
7 Cs 

oj + Cgfc oj — Cgk 

(w + Csfc)2 + (rsfc2)2 _|_ (rsfc2)2 


Cy{k,t) 


{v'^{k,t) ■ v^{-k,0)) ^ 

(t)^(fc, 0 ) 2 ) 


(17) 


where v = y/g is the kinematic viscosity and g the mass 
density. 

Figure [l] shows examples of autocorrelation functions 
for various collision-time steps extracted from simula¬ 
tions without a thermostat. Evidently, the obtained 
Cy(t) decay exponentially. A fit with the exponential 
function 0 yields the viscosities listed in Table |lj These 
values are slightly larger than the theoretical values cal¬ 
culated according to Eqs. ([^ and Q, which is a conse¬ 
quence of the applied approximations in the derivation 
of the theoretical expressions. 

We performed various simulations applying a thermo¬ 
stat and calculated Cy{k,t). Within the accuracy of the 
results, we did not detect any difference between simula¬ 
tions with and without thermostat. 

Eor the sake of completeness, we like to mention 
that Eourier transformation of the correlation function 
Cy{k,t) in Eq. ( [T7| ) yields the well-known long-time tail, 
characteristic for hydrodynamic correlations [571 [Ml 15]. 
Eurther details are presented in Ref. [57]. 


B. Dynamic Structure Factor 


The dynamic structure factor is defined as 
1 /■°° 

S{k,u}) = — {Sp{k,t)dp{-k,0))e (18) 

* J —OO 

in terms of the (number) density fluctuations dp{r,t) = 
p(r, t) — p [m [7^ ini HSl H] , where p denotes the mean 
density and 

N 

p(fc,t)=^e*'=''‘W. (19) 

i=l 

Explicitly, the normalized dynamic structure factor S 
(f Sduj = 1) of an adiabatic system for small k values 


where Cs = a/ yfesE/m is the adiabatic velocity of sound, 
Dt the thermal diffusion coefficient, Fg the sound attenu¬ 
ation factor, and 7 the adiabatic index. More definitions 
and the relation with the MFC parameters are provided 
in Appendix]^ The expression for an isothermal system 
follows by setting Dt = 0 and 7=1 [55] : 

Vk 

2'KS{k,uj) = — 
c 

2ck — OJ 

^ (UJ - ck f + (rfc 2)2 


2ck + OJ 


(oj + ckf + (rfc2)2 


( 21 ) 


Here, c = ^/ksT^m denotes the isothermal speed of 
sound and F the isothermal sound attenuation factor (see 
Appendix 0 . Note that the structure factor is related 
to the longitudinal velocity autocorrelation function via 

[Tiin] 

^ J{v^(k,t)v^(-k,0))e-^‘^^dt= S(k,oj). 

^ ( 22 ) 


This correlation function lacks a Rayleigh line due to the 
appearing frequency (w^) on the right-hand side. 

Figures and 1^ provide examples of S(k, oj) for the 
collision times h — O.lr and 3.Or, respectively, and the 
MBS and MCS scaling schemes. For LSS, we obtain the 
identical structure factors as for MBS within the accu¬ 
racy of the simulations. For the short collision-time step 
(Fig.§, two Brillouin lines are present at the frequen¬ 
cies OJ K, ±cfc. No central Rayleigh line is present, hence, 
there are no temperature fluctuations. The simulation 
result of the MBS thermostat is in very close agreement 
with the theoretical prediction, whereas the height of the 
Brillouin peaks is smaller for the MCS thermostat, but 
the peak positions correspond to those of an isothermal 
system. 

Similarly, for the simulations with h = 3.Or (Fig. |^, 
the Brillouin lines of the MBS thermostat correspond to 
those of an isothermal system, although the peak height 
is somewhat smaller than that of an isothermal system. 
In contrast, the structure factor for the MCS thermo¬ 
stat is close to the theoretical expression of an adibatic 
system. The Brillouin peaks shift to the frequencies 
OJ « zLyGyck, corresponding to adibatic sound propaga¬ 
tion. More importantly, there is a central Rayleigh line. 
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Thus, the MCS thermostat at large collision-time steps 
is not reproducing an isothermal but rather an adiabatic 
system. This may not necessarily be a problem for tem¬ 
perature control, since the Monte Carlo procedure ap¬ 
proaches the desired canonical velocity distribution in 
the limit of a large number of attempts; however, the 
temperature fluctuations are not correct locally. In ad¬ 
dition, typically collision-time steps h < 0.2r are used 
to simulate fluids. Here, a nearly isothermal system is 
achieved for MCS. 

The deviation between the theoretical structure fac¬ 
tor of an isothermal system and the simulation data for 
MCS at h = O.It and MBS at h = 3.Or, respectively, 
indicates that neither method controls temperature per¬ 
fectly locally for these time steps. We attribute the de¬ 
viation from the isothermal dynamic structure factor to 
streaming of the MFC particles. For the MBS thermostat 
and the collision-time step h = O.lr, there is very little 
energy transport during streaming, and thus, the sys¬ 
tem closely resembles an isothermal ensemble. However, 
for h = 3.Or, there is a considerable energy transfer to 
nearby collision cells in the streaming step, which implies 
non-isothermal fluctuations. In case of the MCS method, 
velocity scaling occurs with a certain probability only, 
which leads to large displacements of particles without 
real temperature control. This is particularly pronounced 
for h = 3.Or, where heat is transferred over large dis¬ 
tances during streaming and gives rise to adiabatic rather 
than isothermal fluctuations. The crossover from isother¬ 
mal to adiabatic density fluctuations has been addressed 
in Ref. [55] . 

We finally would like to emphasize that the dynamic 
structure factor for systems with the MCS thermostat 
depends on the parameter Simulations with the “ex¬ 
treme” values ( = 0.05 and 0.3 lead to slight shifts of the 
Brillouin lines and variation in the peak heights. How¬ 
ever, for h = 3t, there is always a pronounced Rayleigh 
peak. 


V. RESULTS: NON-EQUILIBRIUM 
SIMULATIONS 

We determine the fluid viscosity via non-equilibrium 
simulations in order to demonstrate that the viscosity is 
independent of the thermostat and, moreover, that the 
thermostat is not interfering with the flow. 



FIG. 2. (Color online) Normalized dynamic structure factors 
for h = O.lr. The solid line with the smallest peaks (blue) cor¬ 
respond to the MCS thermostat, the line with the next larger 
peaks (red) to the MBS thermostat, and the curve with the 
most pronounced peaks (black) is the theoretical structure 
factor of an isothermal system (Eq. (211). The dashed curve 
indicates the theoretical structure factor of an adiabatic sys¬ 
tem (Eqs. @)- The system size is L/a = 30. 


TABLE I. Kinematic viscosities v = rj/g and their devia¬ 
tions A.V = {v — X 100% from theoretical values 

, Eqs. ^ andj^] obtained from the veloc¬ 
ity autocorrelation function ( |17[ ) and shear-flow simulations 
for the thermalization methods: local simple scaling (LSS), 
Monte Carlo scaling (MCS) [6], and Maxwell-Boltzmann scal¬ 
ing (MBS) [75]. The simulation box size is set to L = 60a in 
the calculation of the VACF for h/r — 1, 2, and 3. 



h/r 

0.1 

0.2 

0.5 

1.0 

2.0 

3.0 

Theory 


0.870 

0.508 

0.407 

0.568 

1.014 

1.486 

VACF 

ulia^lr) 

0.873 

0.515 

0.409 

0.569 

1.006 

1.484 

Aul% 

0.4 

1.3 

0.5 

0.2 

-0.8 

-0.1 

LSS 

i//(a7r) 

0.872 

0.517 

0.411 

0.571 

1.017 

1.492 

Au/% 

0.2 

1.7 

0.9 

0.4 

0.4 

0.4 

MCS 

u/ia^lr) 

0.869 

0.515 

0.412 

0.571 

1.016 

1.493 

Aul% 

-0.1 

1.4 

1.2 

0.4 

0.3 

0.5 

MBS 


0.871 

0.517 

0.414 

0.573 

1.019 

1.494 

Aul% 

0.1 

1.8 

1.5 

0.8 

0.5 

0.5 


A. Shear Elow 

We perform shear-flow simulations for various collision 
times and the LSS, MCS, and MBS thermostat. In all 
cases, we obtain a linear velocity profile, in agreement 
with the theoretical expectation. From the stress ten¬ 
sor values, we calculate the viscosities listed in Table |l| 
The viscosities attained by the various thermostats are 
in close agreement with each other and are in remarkable 


agreement with the theoretical prediction. As expected, 
the simulation values are typically slightly larger than 
the theoretically determined viscosities. However, they 
agree within about 2%. The largest deviation appears for 
h = 0.2r. This supports our expectation that the theo¬ 
retically derived expression for the viscosity agrees well 
with simulation results for larger and smaller collision¬ 
time steps. The agreement between the viscosities ex¬ 
tracted from simulations is even better; the relative error 
is below 1%. 

























Q)/[kgT/(ma^)]^'^ 


FIG. 3. (Color online) Normalized dynamic structure factors 
for h — 3.Or. The solid line with the smallest peaks (blue) cor¬ 
respond to the MCS thermostat, the line with the next larger 
peaks (red) to the MBS thermostat, and the curve with the 
most pronounced peaks (black) is the theoretical structure 
factor of an isothermal system (Eq. (211). The dashed curve 
indicates the theoretical structure factor of an adiabatic sys¬ 
tem (Eqs. (|20[| The system size is L/a = 30. 
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FIG. 4. (Color online) Poiseuille-flow velocity profiles of sim¬ 
ulations with non-zero average phantom-particle momenta, 
i.e., no-slip boundary conditions, for the simple scaling (LSS) 
(short-dashed, red), the Monte Carlo scaling (MCS) (dotted, 
green), and the Maxwell-Boltzmann (MBS) (dashed, blue) 
thermostat. The parabolic velocity profiles with the theoreti¬ 
cally determined viscosities (cf. Table |I| and zero slip length 
are shown by solid lines (black). The collision-time steps are 
h/r = 1.0, 0.1, and 3.0 (top to bottom). 


B. Poiseuille Flow 


Figure shows velocity profiles for a force-driven MFC 
fluid confined between hard walls thermalized by the LSS, 
MCS, or MBS method, and various collision-time steps. 
Here, the average momentum of a phantom particle in a 
wall collision cell is determined by the desired parabolic 
velocity profile (cf. Sec. IID). As displayed in the figure, 


this choice yields a zero fluid velocity at the surface (see 
als Fig. [^. Independent of the applied thermostat, the 
fluid particle temperature and density across the channel 
are constant. For every collision-time step, we find good 
agreement between the velocity profiles of the various 
thermostats. Moreover, the profiles agree well with the 
parabola with the theoretically determined viscosities. 
The actually determined viscosities and their deviations 
from the theoretical values are summarized in Table El 
Here, the simulation data are fitted by the parabola 


Vxiy) = ^{y + ls){L + ls- y), 


(23) 


which yields the slip length 1^ and the kinematic viscosity 
v. As expected, we find a zero slip length in simulations 
where the phantom-particle momentum (P) ^ 0. There 
are only very minor differences between the viscosities 
obtained for the various thermostats, and the viscosi¬ 
ties themselves agree well with the theoretical values. 
Thereby, the numerical values are typically somewhat 
larger, up to approximately 2%. The shear-flow simu¬ 
lations show the same trend. 


Figure compares velocity profiles for no-slip bound¬ 
ary conditions, where (P) 7 ^ 0 , with results with residual 
slip, where we set (P) = 0 (cf. Sec. HD). There is a 
finite slip for (P) = 0, which implies a shift of the whole 
velocity profile to larger velocities. A fit by the profile 
(23) yields the slip length Ig and the viscosity Vg- The 
respective values are summarized in Table |TT1 I n general, 
the profiles are excellently fitted by Eq. (23). Despite 
the differences of the profiles, the viscosities are in close 
agreement. Since a finite residual slip does not alter the 
viscosity, a fit with a finite slip length yields a very ac¬ 
curate estimation of the viscosity. However, a fit with 
zero slip length provides a somewhat different viscosity. 
Thereby, the overall numerical profile is not very well re¬ 
produced by the theoretical parabola. Inclusion of a slip 
length improves fitting considerably. 


VI. DISCUSSION 

The viscosity of the MFC fluid is dominated by con¬ 
tributions from collisions at small, and by kinetic con¬ 
tributions (streaming) at large collision-time steps. This 
suggests that a random shift of the collision lattice can 
be omitted at large collision-time steps [iniiiHiiiH], and 
partially filled collision cells would not matter anymore. 
However, lack of a random shift causes various ambigui¬ 
ties. Without random shift and phantom particles, there 
are only bounce-back interactions during streaming with 
walls, which does not prevent slip strictly, because the 
average velocity at the wall will never be zero during a 
MFC dynamics step. More severely, the induced veloc- 
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TABLE 11. Kinematic viscosities v and their deviations Ai/ 
with respect to theoretical values extracted from Poiseuille 
flow simulations by fitting a parabola to the flow profile for 
the various thermostats, la is the slip length and Va the kine¬ 
matic viscosity for simulations, where the average phantom- 
particle momentum is set to zero. The other viscosities 
are obtain form simulations with a profile-matched phantom- 
particle momentum. 



hjr 

0.1 

1.0 

3.0 


vHa^/r) 

0.886 

0.570 

1.483 


Avl% 

1.9 

0.4 

-0.2 

LSS 

Valia? It) 

0.888 

0.572 

1.486 


la la 

0.167 

0.054 

0.079 


AVal% 

2.0 

0.6 

0.0 


vlia^lr) 

0.869 

0.570 

1.483 


Avl% 

-0.1 

0.2 

-0.2 

MCS 

Vs Ha? It) 

0.869 

0.571 

1.482 


la! a 

0.170 

0.057 

0.093 


AVal% 

-0.1 

0.4 

-0.3 


vHa^h) 

0.882 

0.573 

1.485 


Avl% 

1.4 

0.7 

-0.1 

MBS 

Vs Ha? It) 

0.882 

0.574 

1.488 


Is !Oi 

0.172 

0.057 

0.090 


AVal% 

1.4 

0.9 

0.1 



FIG. 5. (Color online) Velocity profiles of the Poiseuille flow 
with residual slip (upper curve) and without slip (bottom 
curve) by assigning a finite (negative) velocity to phantom 
particles. The time step is h = O.lr. The inset shows the 
profiles close to the solid wall at y = 0. The dashed line indi¬ 
cates the fit of the velocity profile (231 with finite slip length. 



FIG. 6. (Color online) Velocity profiles for the Monte 
Carlo scaling (MCS) approach (top, green) and the Maxwell- 
Boltzmann scaling (MBS) method (bottom, blue) for the 
collision-time step h = 3.Or. The result for the LSS ap¬ 
proach is indistinguishable from the MBS result. No random 
shift is applied. The smooth solid line (black) is the theoreti¬ 
cal parabolic velocity profile with the analytically determined 
viscosity presented in Table [I] 


schemes with conserved linear momentum, since the sta¬ 
tionary state distribution of the relative velocities Avi 
is Gaussian with zero mean. Hence, the average veloc¬ 
ity of a particle in a cell after collision is (vi) = {vcm}- 
The calculation of the velocities after streaming yields 
a smoother profile, in particular for very large collision¬ 
time steps. If only the velocity of the cell center would 
be considered, i.e., the bin width for the calculation of 
the profile is set equal to the size of the collision cell, 
the steps are invisible and a smooth profile is obtained. 
As revealed by the in-depth studies of Ref. [7B], lack or 
presence of a random shift leads to slightly different ve¬ 
locity profiles, with a higher viscosity in the presence of 
a random shift. The difference for the studied time step, 
however, is extremely small and is expected to be even 
smaller for larger h. Importantly, the difference between 
the velocity profiles is not related to partially filled colli¬ 
sion cells, but only to the shift of the collision lattice. To 
avoid ambiguities in the calculation of the velocity profile 
(after streaming versus after collision), we recommend to 
use a random shift of the collision lattice for any time 
step. This yields a unique viscosity. 


VII. SUMMARY AND CONCLUSIONS 


ity profile is no longer smooth on the length scale of a 
collision cell. As shown in Fig. correlations on the 
cell level lead to essentially constant average velocities of 
the particles in a collision cell and, hence, to a step-like 
overall profile. Note that we calculate the velocity pro¬ 
file after a collision. The steps appear for all collision 


We have presented a detailed evaluation of various 
thermostating approaches for non-equilibrium multipar¬ 
ticle collision dynamics simulations. The purpose of our 
studies is twofold. On the one hand, we want to shed light 
on the accuracy with which the non-equilibrium aspects 
of the fluid are reproduced or are perturbed by a partic- 
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ular thermostat. On the other hand, we intend to clarify 
whether there are deviations, and if so, how large, be¬ 
tween the analytically determined viscosity of the MPC 
fluid and that extracted from simulation. For this pur¬ 
pose, we have determined the fluid viscosities for various 
MPC collision-time steps by equilibrium simulations via 
the transverse fluid velocity autocorrelation function and 
by non-equilibrium simulations calculating the stress ten¬ 
sor under simple shear, as well as the velocity proHle for 
a Poiseuille flow. In the non-equilibrium simulations, we 
control temperature on the cell level by three different 
methods: local simple scaling (LSS), a Monte Carlo-like 
scheme (MCS) [7H] . and by the Maxwell-Boltzmann scal¬ 
ing (MBS) approach [75] . 

The calculation of the dynamic structure factor with 
the equilibrium fluid density fluctuations of the system 
thermalized by the MCS method yields, at hrst unex¬ 
pected, a pronounced Rayleigh peak for collision-time 
steps h/r > 1. Hence, in such a system thermal fluc¬ 
tuations play a major role and it is closer to an isoen- 
tropic rather than an isothermal ensemble. The correct 
average asymptotic temperature is assumed for many 
Monte Carlo steps, however, the fluctuations do not cor¬ 
respond to an isothermal ensemble. This implies different 
transport coefficients compared to an strict isothermal 
system—they may neither correspond to an adiabatic nor 
to an isothermal system. 

However, our simulations suggest that every employed 
thermostat leaves the viscosity unchanged, or at least 
affects it to such small extent that it is difficult to de¬ 
tect by the flow profiles or thermodynamic properties. 
Hence, we find only minor differences between the vis¬ 
cosities obtained by the various approaches. Thus, we 
consider all of them suitable for non-equilibrium simula¬ 
tions for weakly perturbed systems. The drawback of the 
LSS method, however, is that the velocity distribution of 
the fluid particles is non-Gaussian, which leads to arti¬ 
facts in the density and even temperature distribution at 
larger field strengths. As discussed in Ref. [75], the MBS 
method provides accurate results even at large fields. 

Looking at the agreement between the analytically 
predicted viscosities with those determined by simula¬ 
tions, we hnd slightly larger numerical values in the range 
0.2 < hjr < 1 than theoretically predicted (cf. Tablejl]). 
All thermostats yield consistently slightly larger viscosi¬ 
ties, with some small variations. 

We have only considered the SRD variant of MPC, 
where fluid velocities are rotated around a randomly ori¬ 
ented axis. In Ref. [7S], other collision rules have been 
applied for Poiseuille flow simulations, in particular ro¬ 
tations around the Cartesian axes only. Simulations 
exploiting the MPC-AT method yield velocity profiles, 
which agree very well with those determined with the the¬ 
oretical viscosity for large collision-time steps. However, 
simulations applying rotations of the relative velocities 
around one of the randomly selected Cartesian axis yields 
considerable deviations between simulation and theoret¬ 
ical results. Applying the same rule, we also find larger 


deviations than those found by the above applied collision 
rule. Hence, the collision rule affects the fluid behavior 
under non-equilibrium conditions. In the axis-rotation 
scheme, there seem to be considerable correlations of the 
fluid particles in a collision cell, more than in the other 
algorithms. 

Our simulations reveal a major effect of the violation 
of Galilean invariance on the flow properties in the form 
of stair-like prohles, for both, simple shear and Poiseuille 
flow. The effect as such is independent of the collision¬ 
time step. Signatures of such steps have also been re¬ 
ported in Ref. m- As we have shown, the steps com¬ 
pletely disappear when a random shift of the collision 
lattice is applied. Thus, we strongly recommend to ap¬ 
ply a random shift of the collision lattice even for large 
collision-time steps, although physically relevant fluid 
properties can only be expected on length scales larger 
than a collision cell. 

Moreover, the random shift is intimately connected 
with the boundary condition. A no-slip boundary condi¬ 
tion is best fulfilled by applying a random shift and in¬ 
clusion of phantom particles nniiiH], for both, stationary 
surfaces as well as dissolved solid bodies liinz]- Further 
investigations of the boundary conditions on the dynam¬ 
ics of colloids are currently under way, with emphasis 
on the differences in colloid dynamics between slip and 
no-slip boundaries. 


Appendix A: Transport Coefficients 


Here, the various transport coefficients defined in Sec¬ 
tion |IV| are given in terms of the MPC-SRD fluid pa¬ 
rameters naiHS]. We assume that (Nc) 1, such that 
+ « (iVc)-l. 

The thermal diffusion coefficient Dt is given by 

Dt = D^t + (A1) 


with 


ibh (w) 


= 


ksTh 

2m 


1 - 
3 


1 


1 — cos a 


(Nc) 
-1-b 


(1 — cos a) , 


(A2) 


(Ne) 


1 1 


4 sin^ a/2 


(A3) 


The specific heat capacities are 


_ 3A:b 
2m 


Cp 




5 

3 


(A4) 


The sound attenuation factor of an adiabatic system is 
defined as 

= + (A5) 

with 

+ = (A6) 
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and the kinematic viscosity v = yy/p = + v'^. The is 

viscosities rj'^ and are defined in Eqs. (|^ and §■ 1 

r = (A7) 

For an isothermal fiuid, the sound attenuation factor 
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